Explore temperature data

Load libraries

Code

pkgs <- c("here", "tidyverse", "readxl", "RColorBrewer", "ISOweek", "viridis", "ncdf4", "reshape2") 

# minpack.lm needed if using nlsLM()
if(length(setdiff(pkgs,rownames(installed.packages()))) > 0){

    install.packages(setdiff(pkgs,rownames(installed.packages())),dependencies = T)
  
  }

invisible(lapply(pkgs, library, character.only = T))

# devtools::install_github("seananderson/ggsidekick") # not on CRAN 
library(ggsidekick)
theme_set(theme_sleek())

# Load functions
home <- here::here()
fxn <- list.files(paste0(home, "/R/functions"))
invisible(sapply(FUN = source, paste0(home, "/R/functions/", fxn)))

Temperature logger data

Code

home <- here::here()
filenames <- list.files(paste0(home,"/data/temp_data"))
raw_data <- NULL

for(i in 1:length(filenames)) {
  dat <- read_excel(paste0(home, "/data/temp_data/", filenames[i]))
  if(i == 1) names <- colnames(dat)
  if(i != 1) names(dat) <- names
  area <- substr(filenames[i], 1, 2)
  if(area == "SI") area <- substr(filenames[i], 1, 5)
  if(area == "JM") area <- paste0(substr(filenames[i], 1, 2), substr(filenames[i], 9, 10))
  if(area == "JMT0") area <- substr(filenames[i], 1, 2)
  dat$area <- area
  raw_data <- data.frame(rbind(raw_data, dat))
} # JMT0 is really JM_T10 and taken as the main temperature time series for JM
#> New names:
#> • `` -> `...3`

# Filter only one JM series: Which one?
raw_data <- raw_data %>% filter(!area %in% c("JMT2", "JMT3"))

# Depth measurements
log_dat <- raw_data %>%
  mutate(Depth = gsub(",", ".", gsub('m', '', Depth))) %>%
  filter(Depth != "Yta") %>%
  filter(Depth != "Siktdjup") %>%
  mutate_at('Depth', ~str_replace(., "0.3-0.5", "0.5")) %>%
  mutate_at('Depth', ~str_replace(., "ca 1.5", "1.5")) %>%
  mutate(Depth = as.numeric(Depth)) %>%
  drop_na(Depth)

# Filter other columns
log_dat <- log_dat %>%
  select(area, Station_Code, Year, Month, Day, Depth, Mean) %>%
  filter(Depth >= 0.5 & Depth <= 1.5) %>% # restrict depth range
  filter(Year %in% seq(1963, 2019, 1)) %>% # complete years only
  filter(!is.na(Station_Code)) %>%
  filter(Month %in% seq(12)) %>%
  filter(Day %in% seq(31)) %>%
  mutate(Mean = as.numeric(Mean)) %>%
  select(-Station_Code) %>%
  mutate(date = as.Date(paste(Year, Month, Day, sep = "-")))

# Calculate some date variables needed when joining data sets and fitting models
log_dat <- log_dat %>% 
  mutate(yday = yday(date),
         month = month(date)) %>% 
  rename(year = Year,
         temp = Mean) %>% 
  dplyr::select(year, area, temp, yday, month, date) %>% 
  mutate(source = "logger",
         date = as.character(date))

# Plot data over year
log_dat %>%
  ggplot(., aes(x = year, y = temp, color = yday)) +
  geom_point(size = 0.25) +
  scale_color_viridis() +
  theme(plot.title = element_text(size = 15, face = "bold")) +
  theme(axis.text = element_text(size = 12), axis.title = element_text(size = 15)) +
  guides(color = guide_legend(override.aes = list(size = 1))) +
  labs(x = "Year", y = "Temperature") +
  facet_wrap(~area) +
  NULL

Code

# Plot data over yearday
log_dat %>%
  ggplot(.,aes(x = yday, y = temp, color = factor(year))) +
  geom_point(size = 0.25) +
  scale_color_viridis(discrete = TRUE) +
  theme(plot.title = element_text(size = 15, face = "bold")) +
  theme(axis.text = element_text(size = 12), axis.title = element_text(size = 15)) +
  theme(axis.text.x = element_text(angle = 90)) +
  guides(color = guide_legend(override.aes = list(size = 1))) +
  labs(x = "Day-of-the-year", y = "Temperature") +
  # facet_grid(cols=vars(area)) +
  facet_wrap(~area) +
  NULL

Temperature at fishing

Code

kul_dat <- read_excel(paste0(home, "/data/temp_data_fishing/temp_at_fishing.xlsx"), skip = 4) |> 
  mutate(Lokal = ifelse(Lokal == "Simpevarp" & Fångstområde == "Ekö", "SI_EK", Lokal),
         Lokal = ifelse(Lokal == "Simpevarp" & Fångstområde == "Hamnefjärden", "SI_HA", Lokal)) |> 
  mutate(area = recode_factor(Lokal,
                              "Kvädöfjärden" = "JM", 
                              "Muskö" = "MU",
                              "Holmön" = "HO",
                              "Forsmark" = "FM",
                              "Finbo, Åland" = "FB",
                              "Råneå" = "RA",
                              "Torhamn, Karlskrona Ö skärgård" = "TH",
                              "Biotestsjön, Forsmark" = "BT",
                              "Vinö" = "VN")) |> 
  rename(temp = Värde,
         year = År, 
         date = `Vittjnings\ndatum`) |> 
  mutate(db = "KUL",
         yday = yday(date),
         month = month(date),
         id = paste(year, area, sep = "_")) |> 
  mutate(temp = ifelse(area == "SI_HA" & temp == 267, 26.7, temp)) |> 
  filter(!area == "Simpevarp") |> 
  filter(temp < 100) |> 
  filter(Vertikalstyp == "Yta") |> 
  select(year, yday, date, temp, area, id, db)
#> New names:
#> • `` -> `...8`
#> • `` -> `...10`

# Plot
ggplot(kul_dat, aes(x = year, y = temp)) +
  geom_point(size = 1) +
  guides(color = "none") +
  labs(x = "Year", y = "Temperature") +
  facet_wrap(~area) +
  NULL

Code

# Now we need to add in the rest of the data from Firre database
home <- here::here()
filenames <- list.files(paste0(home,"/data/temp_data_fishing_firre"))
raw_data<-NULL
for(i in 1:length(filenames)) {
  dat <- read_excel(paste0(home,"/data/temp_data_fishing_firre/",filenames[i]))
  if(i == 1) names <- colnames(dat)
  if(i != 1) names(dat) <- names
  area <- substr(filenames[i], 1, 2)
  if(area == "SI") area <- substr(filenames[i], 1, 5)
  #if(area=="JM") area <- paste0(substr(filenames[i],1,2),substr(filenames[i],9,10))
  if(area == "JMT0") area <- substr(filenames[i], 1, 2)
  dat$area <- area
  raw_data <- data.frame(rbind(raw_data, dat))
} # JMT0 is really JM_T10 and taken as the main temperature time series for JM

# Here I get a warning when I convert to data from Year, Week, day-of-week. Could be due to different date standards (week starting with 0 instead of 1, which I think is the ISO way and how our data are coded. It's however not an option in as.Date, which is why i looked into `ISOweek2date`. But that throws an error for some of the dates. Hmm. A simple test however shows that as.Date returns the correct day (double checked). The warnings from as.Date and the error from ISOweek2date could be due to error in data entry. Will explore. 
# https://stackoverflow.com/questions/45549449/transform-year-week-to-date-object
# https://www.r-bloggers.com/2013/08/date-formats-in-r/
firre_dat <- raw_data |> 
  mutate(group = ifelse(VkDag > 100, "2", "1")) |> 
  mutate(VkDag2 = ifelse(group == "1", paste0(0, VkDag), VkDag)) |>  # add zero before strings with length 2...
  separate(VkDag2, sep = 2, into = c("week", 'day'), extra = 'drop', remove = FALSE) |> 
  mutate(week = as.numeric(week), # to get rid of the 0
         day = as.numeric(day)) |> 
  mutate(date = as.Date(paste(Årtal, week, day, sep = "-"), "%Y-%U-%u")) |> 
  # mutate(week2 = paste0("W", week),
  #        weekdate = paste(Årtal, week2, day, sep = "-")) #|> 
  mutate(#date = ISOweek2date(weekdate),
         yday = yday(date),
         month = month(date)) |> 
  rename(year = Årtal,
         stn_nr = Station,
         section_nr = Sektion) |> 
  filter(!if_all(c(MedelFörTemperatur_i, MedelFörTemperatur_u), is.na)) |> 
  mutate(MedelFörTemperatur_i = ifelse(is.na(MedelFörTemperatur_i), MedelFörTemperatur_u, MedelFörTemperatur_i),
         MedelFörTemperatur_u = ifelse(is.na(MedelFörTemperatur_u), MedelFörTemperatur_i, MedelFörTemperatur_u)) |> 
  mutate(temp = (MedelFörTemperatur_i + MedelFörTemperatur_u) / 2,
         db = "FIRRE",
         id = paste(year, area, sep = "_")) |> 
  select(year, yday, month, day, date, VkDag, temp, area, stn_nr, section_nr, db, id) |> 
  drop_na(date) |> 
  mutate(date = as.character(date))
#> Warning: There were 24 warnings in `mutate()`.
#> The first warning was:
#> ℹ In argument: `date = as.Date(paste(Årtal, week, day, sep = "-"),
#>   "%Y-%U-%u")`.
#> Caused by warning in `strptime()`:
#> ! (0-based) yday 369 in year 1972 is invalid
#> ℹ Run `dplyr::last_dplyr_warnings()` to see the 23 remaining warnings.

# Now we need to merge it with fish_dat
fish_dat <- bind_rows(firre_dat |> filter(!id %in% kul_dat$id), # remove duplicates
                      kul_dat) |> 
  mutate(source = "fishing")

head(fish_dat)
#>   year yday month day       date VkDag  temp area stn_nr section_nr    db
#> 1 1991  239     8   2 1991-08-27   342 15.75   BS      1          1 FIRRE
#> 2 1992  224     8   2 1992-08-11   322 16.65   BS      1          1 FIRRE
#> 3 1992  226     8   4 1992-08-13   324 16.60   BS      1          1 FIRRE
#> 4 1992  228     8   6 1992-08-15   326 15.85   BS      1          1 FIRRE
#> 5 1994  266     9   5 1994-09-23   385 12.30   BS      1          1 FIRRE
#> 6 1991  239     8   2 1991-08-27   342 15.65   BS      2          1 FIRRE
#>        id  source
#> 1 1991_BS fishing
#> 2 1992_BS fishing
#> 3 1992_BS fishing
#> 4 1992_BS fishing
#> 5 1994_BS fishing
#> 6 1991_BS fishing

fish_dat |> 
  filter(area == "BT") |> 
  group_by(source) |> 
  summarise(min_yr = min(year))
#> # A tibble: 1 × 2
#>   source  min_yr
#>   <chr>    <dbl>
#> 1 fishing   1984

# Plot data over year
fish_dat |>
  ggplot(aes(x = year, y = temp, color = yday)) +
  geom_point(size = 0.25) +
  scale_color_viridis() +
  theme(plot.title = element_text(size = 15, face = "bold")) +
  theme(axis.text.x = element_text(angle = 90)) +
  theme(axis.text = element_text(size = 12), axis.title = element_text(size = 15)) +
  guides(color = guide_legend(override.aes = list(size = 1))) +
  labs(x = "Year", y = "Temperature") +
  facet_wrap(~area) +
  NULL

Code

# Plot data over yearday
fish_dat |>
  ggplot(aes(x = yday, y = temp, color = factor(year))) +
  geom_point(size = 0.25) +
  scale_color_viridis(discrete = TRUE) +
  theme(plot.title = element_text(size = 15, face = "bold")) +
  theme(axis.text = element_text(size = 12), axis.title = element_text(size = 15)) +
  theme(axis.text.x = element_text(angle = 90)) +
  guides(color = guide_legend(override.aes = list(size = 1))) +
  labs(x = "Day-of-the-year", y = "Temperature") +
  facet_wrap(~area) +
  NULL

ERSST data

Code

# Longitude and latitude attributes for each area
area <- c("BS", "BT", "FB", "FM", "HO", "JM", "MU", "RA", "SI_EK", "SI_HA", "TH", "VN")
nareas <- length(area)
lat <- c(60, 60.4, 60.3, 60.5, 63.7, 58, 59, 65.9, 57.3, 57.4, 56.1, 57.5)
lon <- c(21.5, 18.1, 19.5, 18, 20.9, 16.8, 18.1, 22.3, 16.6, 16.7, 15.9, 16.9)
area_attr <- data.frame(cbind(area = area, lat = lat, lon = lon)) |>
  mutate_at(c("lat", "lon"), as.numeric)

# Download ERSST data - only run when update needed
# sst_dat_all <- get_ersst_data(years=seq(1940,2022),data.dir=paste0(home,"/data"), latrange=c(55,66),lonrange=c(15,23))

# SST based on ERSST data with relatively low spatial resolution (2x2 degrees)
# Need to cover at least 2x2 grid area with even numbers for longitude/latitude
sst_areas <- NULL
lat_ranges <- lon_ranges <- list() # for testing only
for(a in 1:nareas) {
  lat_range <- c(2*floor(area_attr$lat[a]/2), 2*floor(area_attr$lat[a]/2) + 2)
  lon_range <- c(2*floor(area_attr$lon[a]/2), 2*floor(area_attr$lon[a]/2) + 2)
  sst_area <- load_ersst_data(years = c(1940, 2022), data.dir = paste0(home, "/data"), ncfilename = "sst.mnmean.nc", latrange = lat_range, lonrange = lon_range)
  sst_area$area <- area_attr$area[a]
  sst_areas <- bind_rows(sst_areas, sst_area)
  lat_ranges[[a]] <- lat_range
  lon_ranges[[a]] <- lon_range
}
#> Joining with `by = join_by(lon_index)`
#> Joining with `by = join_by(lat_index)`
#> Joining with `by = join_by(time_index)`
#> Rows: 3944 Columns: 7
#> ── Column specification
#> ──────────────────────────────────────────────────────── Delimiter: "," chr
#> (1): month dbl (5): sst, lon, lat, time, year date (1): date
#> ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
#> Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> `summarise()` has grouped output by 'year'. You can override using the
#> `.groups` argument.
#> Joining with `by = join_by(lon_index)`
#> Joining with `by = join_by(lat_index)`
#> Joining with `by = join_by(time_index)`
#> Rows: 3944 Columns: 7
#> ── Column specification
#> ──────────────────────────────────────────────────────── Delimiter: "," chr
#> (1): month dbl (5): sst, lon, lat, time, year date (1): date
#> ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
#> Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> `summarise()` has grouped output by 'year'. You can override using the
#> `.groups` argument.
#> Joining with `by = join_by(lon_index)`
#> Joining with `by = join_by(lat_index)`
#> Joining with `by = join_by(time_index)`
#> Rows: 3944 Columns: 7
#> ── Column specification
#> ──────────────────────────────────────────────────────── Delimiter: "," chr
#> (1): month dbl (5): sst, lon, lat, time, year date (1): date
#> ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
#> Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> `summarise()` has grouped output by 'year'. You can override using the
#> `.groups` argument.
#> Joining with `by = join_by(lon_index)`
#> Joining with `by = join_by(lat_index)`
#> Joining with `by = join_by(time_index)`
#> Rows: 3944 Columns: 7
#> ── Column specification
#> ──────────────────────────────────────────────────────── Delimiter: "," chr
#> (1): month dbl (5): sst, lon, lat, time, year date (1): date
#> ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
#> Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> `summarise()` has grouped output by 'year'. You can override using the
#> `.groups` argument.
#> Joining with `by = join_by(lon_index)`
#> Joining with `by = join_by(lat_index)`
#> Joining with `by = join_by(time_index)`
#> Rows: 3944 Columns: 7
#> ── Column specification
#> ──────────────────────────────────────────────────────── Delimiter: "," chr
#> (1): month dbl (5): sst, lon, lat, time, year date (1): date
#> ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
#> Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> `summarise()` has grouped output by 'year'. You can override using the
#> `.groups` argument.
#> Joining with `by = join_by(lon_index)`
#> Joining with `by = join_by(lat_index)`
#> Joining with `by = join_by(time_index)`
#> Rows: 3944 Columns: 7
#> ── Column specification
#> ──────────────────────────────────────────────────────── Delimiter: "," chr
#> (1): month dbl (5): sst, lon, lat, time, year date (1): date
#> ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
#> Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> `summarise()` has grouped output by 'year'. You can override using the
#> `.groups` argument.
#> Joining with `by = join_by(lon_index)`
#> Joining with `by = join_by(lat_index)`
#> Joining with `by = join_by(time_index)`
#> Rows: 3944 Columns: 7
#> ── Column specification
#> ──────────────────────────────────────────────────────── Delimiter: "," chr
#> (1): month dbl (5): sst, lon, lat, time, year date (1): date
#> ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
#> Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> `summarise()` has grouped output by 'year'. You can override using the
#> `.groups` argument.
#> Joining with `by = join_by(lon_index)`
#> Joining with `by = join_by(lat_index)`
#> Joining with `by = join_by(time_index)`
#> Rows: 3944 Columns: 7
#> ── Column specification
#> ──────────────────────────────────────────────────────── Delimiter: "," chr
#> (1): month dbl (5): sst, lon, lat, time, year date (1): date
#> ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
#> Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> `summarise()` has grouped output by 'year'. You can override using the
#> `.groups` argument.
#> Joining with `by = join_by(lon_index)`
#> Joining with `by = join_by(lat_index)`
#> Joining with `by = join_by(time_index)`
#> Rows: 3944 Columns: 7
#> ── Column specification
#> ──────────────────────────────────────────────────────── Delimiter: "," chr
#> (1): month dbl (5): sst, lon, lat, time, year date (1): date
#> ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
#> Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> `summarise()` has grouped output by 'year'. You can override using the
#> `.groups` argument.
#> Joining with `by = join_by(lon_index)`
#> Joining with `by = join_by(lat_index)`
#> Joining with `by = join_by(time_index)`
#> Rows: 3944 Columns: 7
#> ── Column specification
#> ──────────────────────────────────────────────────────── Delimiter: "," chr
#> (1): month dbl (5): sst, lon, lat, time, year date (1): date
#> ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
#> Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> `summarise()` has grouped output by 'year'. You can override using the
#> `.groups` argument.
#> Joining with `by = join_by(lon_index)`
#> Joining with `by = join_by(lat_index)`
#> Joining with `by = join_by(time_index)`
#> Rows: 3944 Columns: 7
#> ── Column specification
#> ──────────────────────────────────────────────────────── Delimiter: "," chr
#> (1): month dbl (5): sst, lon, lat, time, year date (1): date
#> ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
#> Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> `summarise()` has grouped output by 'year'. You can override using the
#> `.groups` argument.
#> Joining with `by = join_by(lon_index)`
#> Joining with `by = join_by(lat_index)`
#> Joining with `by = join_by(time_index)`
#> Rows: 3944 Columns: 7
#> ── Column specification
#> ──────────────────────────────────────────────────────── Delimiter: "," chr
#> (1): month dbl (5): sst, lon, lat, time, year date (1): date
#> ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
#> Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> `summarise()` has grouped output by 'year'. You can override using the
#> `.groups` argument.

latranges <- data.frame(matrix(unlist(lat_ranges), ncol = 2, byrow = T))
lonranges <- data.frame(matrix(unlist(lon_ranges), ncol = 2, byrow = T))

# Plot SST by area in each month
sst_areas |>
  ggplot(aes(x = year, y = meanSST, group = as.factor(month), color = as.factor(month))) +
  geom_line() +
  scale_color_viridis(discrete = TRUE) +
  scale_x_continuous(breaks = seq(1940, 2020, 10)) +
  theme(plot.title = element_text(size = 15, face = "bold")) +
  theme(axis.text.x = element_text(angle= 90)) +
  theme(axis.text = element_text(size = 12), axis.title = element_text(size = 15)) +
  labs(x = "Year", y = "Mean SST", title = "Mean SST in each month by area", color = "Month") +
  facet_wrap(~area, scale = "free_y") +
  NULL

Code

sst_areas <- sst_areas |> 
  mutate(date = paste(year, month, 15, sep = "-"),
         yday = yday(date),
         source = "errs",
         month = as.numeric(month)) |> 
  rename(temp = meanSST)

errs_dat <- sst_areas

# Plot SST by month in each area
errs_dat |>
  ggplot(aes(x = year, y = temp, group = as.factor(area), color = as.factor(area))) +
  geom_line() +
  scale_color_brewer(palette = "Set3") +
  scale_x_continuous(breaks = seq(1940, 2020, 10)) +
  theme(plot.title = element_text(size = 15, face = "bold")) +
  theme(axis.text.x = element_text(angle = 90)) +
  theme(axis.text = element_text(size = 12), axis.title = element_text(size = 15)) +
  labs(x = "Year", y = "Mean SST", title = "Mean SST in each month by area", color = "Area") +
  facet_wrap(~month, scale = "free_y") +
  NULL

Code

all_temp <- bind_rows(log_dat, errs_dat, fish_dat)

# See if we can compare the mean temperature in a time of the year that have overlap
all_temp |> 
  filter(area == "FB") |> 
  ggplot(aes(yday, fill = source)) + 
  facet_wrap(~year, scales = "free_y") + 
  geom_density(alpha = 0.5, color = NA)

Code

# Plot some overlapping year and areas
all_temp |> 
  filter(area == "VN" &  year > 1995) |> 
  ggplot(aes(yday, temp, color = source)) +
  geom_point() + 
  facet_wrap(~year)

Code

all_temp |> 
  filter(area == "JM" & year < 1995 & year > 1975) |> 
  ggplot(aes(yday, temp, color = source)) +
  geom_point() + 
  facet_wrap(~year)

Code

all_temp |> 
  filter(area == "FM" & year < 2000 & year > 1980) |> 
  ggplot(aes(yday, temp, color = source)) +
  geom_point() + 
  facet_wrap(~year)

Code

# Some examples: OK match but not in BT which makes sense because the ERSS data isn't fine scale enough
all_temp |> 
  filter(area == "FB" & year %in% c(1998:2001)) |> 
  ggplot(aes(yday, temp, color = source)) +
  geom_point() + 
  geom_line() + 
  facet_wrap(~year)

Code

all_temp |> 
  filter(area == "SI_EK" & year %in% c(1998:2001)) |> 
  ggplot(aes(yday, temp, color = source)) +
  geom_point() + 
  geom_line() + 
  facet_wrap(~year)

Code

all_temp |> 
  filter(area == "BT" & year %in% c(1998:2001)) |> 
  ggplot(aes(yday, temp, color = source)) +
  geom_point() + 
  geom_line() + 
  facet_wrap(~year)

Code

# Take a close look at SI_EK 1995
all_temp |> 
  filter(source == "logger" & area == "SI_EK" & year == 1996) |> 
  ggplot(aes(yday, temp)) +
  geom_point() + 
  facet_grid(area~year)

Code

# Filter the row of very low values - broken logger?
all_temp |> 
  filter(source == "logger" & area == "SI_EK" & year == 1996) |> 
  mutate(keep = ifelse(temp < 6.54, "N", "Y")) |> 
  ggplot(aes(yday, temp, color = keep)) +
  geom_point() + 
  facet_grid(area~year)

Code

all_temp <- all_temp |>
  mutate(keep = ifelse(source == "logger" & area == "SI_EK" & year == 1996 & temp < 6.54, "N", "Y")) |>
  filter(keep == "Y") |>
  select(-keep)

Save data!

Code

# Save data for model fitting
write_csv(all_temp, paste0(home, "/output/temp-data-for-fitting.csv"))